home *** CD-ROM | disk | FTP | other *** search
/ Amiga Format CD 46 / Amiga Format CD46 (1999-10-20)(Future Publishing)(GB)[!][issue 1999-12].iso / -in_the_mag- / reader_requests / scilab / demos / lmitool / phase1.sci < prev    next >
Text File  |  1999-09-16  |  3KB  |  115 lines
  1. function [x,Z,z,ul,iters] = phase1(F,blck_szs,M,nu,abstol,maxiters);
  2. // [x,Z,z,ul,iters] = phase1(F,blck_szs,M,nu,abstol,maxiters);
  3. //
  4. // Find an x s.t. F(x) > 0 and Tr F(x) < M 
  5. //    or prove that no such x exists.
  6.  
  7. // minimize    t 
  8. // subject to  F(x) + t*I = F_0 + x_1*F_1 + ... + x_m*F_m + t*I >= 0 
  9. //             Tr F(x) <= M
  10. //
  11. // maximize    -Tr F_0*(Z-z*I) - M*z
  12. // subject to  0 = Tr F_i*(Z - z*I)
  13. //             1 = Tr Z 
  14. //             Z >= 0, z >= 0
  15. //
  16. // Convergence criterion:
  17. // (1) maxiters is exceeded
  18. // (2) duality gap is less than abstol
  19. // (3) primal objective is less than zero or dual objective is greater
  20. //     than zero
  21. //  
  22. // Input arguments:
  23. // F:         (sum_i n_i^2) x (m+1) matrix
  24. //            [ vec(F_0(1)) vec(F_1(1)) ... vec(F_m(1)) ]
  25. //            [ vec(F_0(2)) vec(F_1(2)) ... vec(F_m(2)) ]
  26. //                ...        ...              ...
  27. //            [ vec(F_0(L)) vec(F_1(L)) ... vec(F_m(L)) ]
  28. //            where F_i(j) denotes the jth block of F_i, 
  29. //            F_i(j) has size n_i x n_i
  30. // blck_szs:  L-vector [n_1 ... n_L], contains dimensions of blocks
  31. // M:         scalar
  32. // nu:        >= 1.0
  33. // abstol:    absolute tolerance
  34. // maxiters:  maximum number of iterations
  35. //
  36. // Output arguments:
  37. // x:         m-vector: 
  38. // Z:         block-diagonal matrix stored as 
  39. //            [ vec(Z(1)); vec(Z(2)); ... ; vec(Z(L)) ]
  40. //            where Z(i) is the ith block
  41. // z:         scalar part of last dual iterate, 
  42. // ul:        2-vector, primal and dual objective
  43. // iters:     scalar
  44. if type(F)=5 then F=full(F);end
  45. [mf,nf]=size(F)
  46. m = nf-1;
  47. if (mf ~= sum(blck_szs^2)) 
  48.    error('Dimensions of F do not match blck_szs.');
  49. end;
  50.  
  51. // mineigF is the smallest eigenvalue of F_0
  52. mineigF = 0.0;
  53. blck_szs=matrix(blck_szs,1,prod(size(blck_szs)));
  54. k=0; for n=blck_szs,
  55.    mineigF = min(mineigF, min(real(spec(matrix(F(k+[1:n*n],1),n,n)))));
  56.    k=k+n*n;   // k = sum n_i*n_i 
  57. end;
  58.  
  59. // I is the identity
  60. I = zeros(mf,1);  
  61. k=0; for n=blck_szs,
  62.    I(k+[1:n*n]) = matrix(eye(n,n),n*n,1);   // identity
  63.    k=k+n*n;   // k = sum n_i*n_i 
  64. end;
  65.  
  66. if (M < I'*F(:,1)+1e-5), 
  67.    error('M must be strictly greater than the trace of F_0.'); 
  68. end;
  69.  
  70. // initial x0 
  71. x0 = [zeros(m,1); max(-1.1*mineigF, 1e-5)];
  72.  
  73. // Z0 is the projection of I on the space Tr F_i*Z = 0
  74. Z0 = I - F(:,2:m+1) * ( F(:,2:m+1) \ I );
  75. // mineigZ is the smallest eigenvalue of Z0
  76. mineigZ = 0.0;
  77. k=0; for n=blck_szs,
  78.    mineigZ = min( mineigZ, min(real(spec(matrix(Z0(k+[1:n*n]),n,n)))) );
  79.    k=k+n*n;   // k = sum n_i*n_i 
  80. end;
  81.  
  82. // z = max ( 1e-5, -1.1*mineigZ)
  83. Z0(k+1) = max( -1.1 * mineigZ, 1e-5 );  // z  
  84. Z0(1:k) = Z0(1:k) + Z0(k+1)*I; 
  85. Z0 = Z0 / (I'*Z0(1:k));    // make Tr Z0 = 1
  86. // add scalar block Tr F(x) <= M
  87.  
  88. F = [F, I; M-I'*F(:,1), -I'*F(:,2:m+1), 0]; 
  89. blck_szs = [blck_szs,1];
  90. c = [zeros(m,1); 1];
  91. //Pack Z0 and F
  92. Z0=pack(Z0,blck_szs);F=pack(F,blck_szs);
  93. pause;
  94. [x,Z,ul,info]=semidef(x0,Z0,F,blck_szs,c,[nu,abstol,-1,0,maxiters]);
  95. nz=prod(size(Z))
  96. z=Z(nz)
  97. Z=unpack(Z(1:nz-1),blck_szs(1:prod(size(blck_szs))-1))
  98. Z = Z(1:k);
  99. x = x(1:m);
  100. iters = info(2);
  101. select info(1)
  102. case 1
  103.      error('Max. iters. exeeded');
  104. case 2 then
  105.      disp('Absolute accuracy reached');
  106. case 3 then
  107.      disp('relative accuracy reached');
  108. case 4 then
  109.      disp('target value reached');
  110. case 5 then
  111.      error('target value not achievable');
  112.        else
  113.      error('semidef fails');
  114. end
  115.